!-------------------------------------------------------------------
! Random number tools and wrappers for the Kruununhaka basis set kit
!
! Copyright (C) 2004-2008 Pekka Manninen
! Copyright (C) 2010 Jussi Lehtola
!
!-------------------------------------------------------------------
! This program is distributed under the terms of the GNU General
! Public License as published by the Free Software Foundation;
! either version 2 of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
!-------------------------------------------------------------------

module random
  use definitions
  implicit none
  private
  ! Get a random integer between [1,n]
  public :: rnd_idx
  ! Get n uniform random numbers on [0,1)
  public :: rnd_vec 
  ! Get a (n,m) matrix of uniform random numbers on [0,1)
  public :: rnd_mat 
  ! Get three random indices for differential evolution algorithm
  public :: rnd_triple
contains

  INTEGER FUNCTION rnd_idx(n)
    ! returns a random integer number between 1,2,...n
    IMPLICIT NONE
    REAL(KIND=prec) rnd_uniform ! function in gslrand.c
    INTEGER, INTENT(IN) :: n
    REAL(KIND=prec) :: x
    x=rnd_uniform()
    rnd_idx = CEILING(REAL(n,prec)*x + 1.0_prec)
  END FUNCTION rnd_idx

  FUNCTION rnd_vec(n)
    !returns a n-dimensional vector of random numbers [0,1[
    IMPLICIT NONE
    REAL(KIND=prec) rnd_uniform ! function in gslrand.c
    INTEGER, INTENT(IN) :: n
    REAL(KIND=prec), DIMENSION(n) :: rnd_vec
    INTEGER :: i

    DO i=1,n
       rnd_vec(i)=rnd_uniform()
    END DO

  END FUNCTION rnd_vec

  FUNCTION rnd_mat(n,m)
    !returns a n,m-dimensional matrix of random numbers [0,1[
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n, m
    REAL(KIND=prec) rnd_uniform ! function in gslrand.c
    REAL(KIND=prec), DIMENSION(n,m) :: rnd_mat
    INTEGER :: i, j
    DO i=1,n
       DO j=1,m
          rnd_mat(i,j)=rnd_uniform()
       END DO
    END DO
  END FUNCTION rnd_mat


  SUBROUTINE rnd_triple(i,n,i1,i2,i3)
    !calculates the indices for {x_a, x_b, x_c}
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, n
    INTEGER, INTENT(OUT) :: i1, i2, i3
    DO
       i1 = rnd_idx(n)
       IF (i1 /= i) EXIT
    END DO
    DO
       i2 = rnd_idx(n)
       IF (i2 /= i .AND. i2 /= i1) EXIT
    END DO
    DO
       i3 = rnd_idx(n)
       IF (ALL(i3 /= (/i,i1,i2/))) EXIT
    END DO
  END SUBROUTINE rnd_triple

end module random
